# Following code provided by Preacher and Selig (2012)

a1 =  0.141 # estimated coefficient JS on grit
a2 =  0.121 # estimated coefficient JS on peer affiliation
a3 =  0.158 # estimated coefficient JS on self-control
a4 =  0.023 # estimated coefficient JS on behavior
a5 =  0.104 # estimated coefficient JS on spiritual
a6 =  0.207 # estimated coefficient JS on self-identity

a1std = 0.0687 # SE of coefficient a1
a2std = 0.0914 # SE of coefficient a2
a3std = 0.0839 # SE of coefficient a3
a4std = 0.0520 # SE of coefficient a4
a5std = 0.0620 # SE of coefficient a5
a6std = 0.0857 # SE of coefficient a6

third.b1 =  0.0570 # estimated coefficient grit on best in third grade
third.b2 =  -0.0252 # estimated coefficient peer affiliation on best in third grade
third.b3 =  0.0257 # estimated coefficient self-control on best in third grade
third.b4 =  -0.0559 # estimated coefficient behavior on best in third grade
third.b5 =  0.0550 # estimated coefficient spiritual on best in third grade
third.b6 =  -0.0238 # estimated coefficient self-identity on best in third grade

third.b1std = 0.0299 # SE on coefficient b1
third.b2std = 0.0300 # SE on coefficient b2
third.b3std = 0.0279 # SE on coefficient b3
third.b4std = 0.0322 # SE on coefficient b4
third.b5std = 0.0364 # SE on coefficient b5
third.b6std = 0.0247 # SE on coefficient b6

elem.b1 = -0.00845 # estimated coefficient grit on best in elementary
elem.b2 = -0.0294 # estimated coefficient peer affiliation on best in elementary
elem.b3 = 0.0742 # estimated coefficient self-control on best in elementary
elem.b4 = -0.00767 # estimated coefficient behavior on best in elementary
elem.b5 = 0.0480 # estimated coefficient spiritual on best in elementary
elem.b6 = 0.0181 # estimated coefficient self-identity on best in elementary

elem.b1std =  0.0282 # SE on coefficient b1
elem.b2std =  0.0327 # SE on coefficient b2
elem.b3std =  0.0243 # SE on coefficient b3
elem.b4std =  0.0377 # SE on coefficient b4
elem.b5std =  0.0428 # SE on coefficient b5
elem.b6std =  0.0254 # SE on coefficient b6

top.b1 =  -0.00115 # estimated coefficient of grit on probability placed in top section
top.b2 =  0.0287 # estimated coefficient of peer-affiliation on probability placed in top section
top.b3 =  0.0306 # estimated coefficient of self-control on probability placed in top section
top.b4 =  -0.000725 # estimated coefficient of behavior on probability placed in top section
top.b5 =  0.0271 # estimated coefficient of spiritual on probability placed in top section
top.b6 =  0.0166 # estimated coefficient of self-identity on probability placed in top section

top.b1std = 0.0220 # SE on coefficient b1
top.b2std = 0.0235 # SE on coefficient b2
top.b3std = 0.0227 # SE on coefficient b3
top.b4std = 0.0278 # SE on coefficient b4
top.b5std = 0.0275 # SE on coefficient b5
top.b6std = 0.0187 # SE on coefficient b6

enroll.b1 = 0.0148 # estimated coefficient of grit on current enrollment
enroll.b2 = 0.0029 # estimated coefficient of peer affiliation on current enrollment  
enroll.b3 = 0.0187 # estimated coefficient of self control on current enrollment
enroll.b4 = -0.0290 # estimated coefficient of behavior on current enrollment
enroll.b5 = 0.0303 # estimated coefficient of self-identity on current enrollment
enroll.b6 = -0.00563 # estimated coefficient of self-idntity on current enrollment

enroll.b1std = 0.0107  # SE on coefficent b1
enroll.b2std = 0.0158  # SE on coefficent b2
enroll.b3std = 0.0164  # SE on coefficent b3
enroll.b4std = 0.0179  # SE on coefficent b4
enroll.b5std = 0.0198  # SE on coefficent b5
enroll.b6std = 0.0116  # SE on coefficent b6

rep = 20000 # number of simulated values
conf = 95 # confidence level

a1vec = rnorm(rep)*a1std+a1 # create vector of simulated a1 coefficients
a2vec = rnorm(rep)*a2std+a2 # create vector of simulated a2 coefficients
a3vec = rnorm(rep)*a3std+a3 # create vector of simulated a3 coefficients
a4vec = rnorm(rep)*a4std+a4 # create vector of simulated a4 coefficients
a5vec = rnorm(rep)*a5std+a5 # create vector of simulated a5 coefficients
a6vec = rnorm(rep)*a6std+a6 # create vector of simulated a6 coefficients

third.b1vec = rnorm(rep)*third.b1std+third.b1 # create vector of simulated b1 coefficients
third.b2vec = rnorm(rep)*third.b2std+third.b2 # create vector of simulated b2 coefficients
third.b3vec = rnorm(rep)*third.b3std+third.b3 # create vector of simulated b3 coefficients 
third.b4vec = rnorm(rep)*third.b4std+third.b4 # create vector of simulated b4 coefficients
third.b5vec = rnorm(rep)*third.b5std+third.b5 # create vector of simulated b5 coefficients
third.b6vec = rnorm(rep)*third.b6std+third.b6 # create vector of simulated b6 coefficients

elem.b1vec = rnorm(rep)*elem.b1std+elem.b1
elem.b2vec = rnorm(rep)*elem.b2std+elem.b2 
elem.b3vec = rnorm(rep)*elem.b3std+elem.b3
elem.b4vec = rnorm(rep)*elem.b4std+elem.b4
elem.b5vec = rnorm(rep)*elem.b5std+elem.b5
elem.b6vec = rnorm(rep)*elem.b6std+elem.b6

top.b1vec = rnorm(rep)*top.b1std+top.b1
top.b2vec = rnorm(rep)*top.b2std+top.b2 
top.b3vec = rnorm(rep)*top.b3std+top.b3
top.b4vec = rnorm(rep)*top.b4std+top.b4
top.b5vec = rnorm(rep)*top.b5std+top.b5
top.b6vec = rnorm(rep)*top.b6std+top.b6

enroll.b1vec = rnorm(rep)*enroll.b1std+enroll.b1
enroll.b2vec = rnorm(rep)*enroll.b2std+enroll.b2 
enroll.b3vec = rnorm(rep)*enroll.b3std+enroll.b3
enroll.b4vec = rnorm(rep)*enroll.b4std+enroll.b4
enroll.b5vec = rnorm(rep)*enroll.b5std+enroll.b5
enroll.b6vec = rnorm(rep)*enroll.b6std+enroll.b6

# Grit on best in third grade
grit_third = (a1vec*third.b1vec)

low = (1-conf/100)/2
upp = ((1-conf/100)/2)+(conf/100)
LL = quantile(grit_third,low) # lower limit of confidence interval
UL = quantile(grit_third,upp) # upper quantile for confidence interval
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
grit.third.fig <- hist(grit_third,breaks = 'FD',col = 'skyblue', 
                           xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                           main = '')

# Peer affiliation on best in third grade
affiliation_third = (a2vec*third.b2vec)

LL = quantile(affiliation_third,low)
UL = quantile(affiliation_third,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
affiliation.third.fig <- hist(affiliation_third,breaks = 'FD',col = 'skyblue',
                                    xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                                    main = '')

# Self-control on best in third grade
selfcontrol_third = (a3vec*third.b3vec)

LL = quantile(selfcontrol_third,low)
UL = quantile(selfcontrol_third,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
selfcontrol.third.fig <- hist(selfcontrol_third,breaks = 'FD',col = 'skyblue',
                              xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                              main = '')

# behavior on best in third grade
behavior_third = (a4vec*third.b4vec)

LL = quantile(behavior_third,low)
UL = quantile(behavior_third,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
behavior.third.fig <- hist(behavior_third,breaks = 'FD',col = 'skyblue',
                              xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                              main = '')

# Spiritual on best in third grade
spiritual_third = (a5vec*third.b5vec)

LL = quantile(spiritual_third,low)
UL = quantile(spiritual_third,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
spiritual.third.fig <- hist(spiritual_third,breaks = 'FD',col = 'skyblue',
                              xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                              main = '')

# Self identity on best in third grade
identity_third = (a6vec*third.b6vec)

LL = quantile(identity_third,low)
UL = quantile(identity_third,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
identity.third.fig <- hist(identity_third,breaks = 'FD',col = 'skyblue',
                              xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                              main = '')

# Grit on best in elementary
grit_elem = (a1vec*elem.b1vec)

LL = quantile(grit_elem,low) # lower limit of confidence interval
UL = quantile(grit_elem,upp) # upper quantile for confidence interval
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
grit.elem.fig <- hist(grit_elem,breaks = 'FD',col = 'skyblue', 
                       xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                       main = '')

# Peer affiliation on best in elementary
affiliation_elem = (a2vec*elem.b2vec)

LL = quantile(affiliation_elem,low)
UL = quantile(affiliation_elem,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
affiliation.elem.fig <- hist(affiliation_elem,breaks = 'FD',col = 'skyblue',
                              xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                              main = '')

# Self-control on best in elementary
selfcontrol_elem = (a3vec*elem.b3vec)

LL = quantile(selfcontrol_elem,low)
UL = quantile(selfcontrol_elem,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
selfcontrol.elem.fig <- hist(selfcontrol_elem,breaks = 'FD',col = 'skyblue',
                              xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                              main = '')

# behavior on best in elementary
behavior_elem = (a4vec*elem.b4vec)

LL = quantile(behavior_elem,low)
UL = quantile(behavior_elem,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
behavior.elem.fig <- hist(behavior_elem,breaks = 'FD',col = 'skyblue',
                           xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                           main = '')

# Spiritual on best in elementary
spiritual_elem = (a5vec*elem.b5vec)

LL = quantile(spiritual_elem,low)
UL = quantile(spiritual_elem,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
spiritual.elem.fig <- hist(spiritual_elem,breaks = 'FD',col = 'skyblue',
                            xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                            main = '')

# Self identity on best in elementary
identity_elem = (a6vec*elem.b6vec)

LL = quantile(identity_elem,low)
UL = quantile(identity_elem,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
identity.elem.fig <- hist(identity_elem,breaks = 'FD',col = 'skyblue',
                           xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                           main = '')


# Grit on probability in top section
grit_top = (a1vec*top.b1vec)

LL = quantile(grit_top,low) # lower limit of confidence interval
UL = quantile(grit_top,upp) # upper quantile for confidence interval
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
grit.top.fig <- hist(grit_top,breaks = 'FD',col = 'skyblue', 
                      xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                      main = '')

# Peer affiliation on probability in top section
affiliation_top = (a2vec*top.b2vec)

LL = quantile(affiliation_top,low)
UL = quantile(affiliation_top,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
affiliation.top.fig <- hist(affiliation_top,breaks = 'FD',col = 'skyblue',
                             xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                             main = '')

# Self-control on probability in top section
selfcontrol_top = (a3vec*top.b3vec)

LL = quantile(selfcontrol_top,low)
UL = quantile(selfcontrol_top,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
selfcontrol.top.fig <- hist(selfcontrol_top,breaks = 'FD',col = 'skyblue',
                             xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                             main = '')

# behavior on probability in top section
behavior_top = (a4vec*top.b4vec)

LL = quantile(behavior_top,low)
UL = quantile(behavior_top,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
behavior.top.fig <- hist(behavior_top,breaks = 'FD',col = 'skyblue',
                          xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                          main = '')

# Spiritual on probability in top section
spiritual_top = (a5vec*top.b5vec)

LL = quantile(spiritual_top,low)
UL = quantile(spiritual_top,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
spiritual.top.fig <- hist(spiritual_top,breaks = 'FD',col = 'skyblue',
                           xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                           main = '')

# Self identity on probability in top section
identity_top = (a6vec*top.b6vec)

LL = quantile(identity_top,low)
UL = quantile(identity_top,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
identity.top.fig <- hist(identity_top,breaks = 'FD',col = 'skyblue',
                          xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                          main = '')



# Grit on currently enrolled
grit_enroll = (a1vec*enroll.b1vec)

LL = quantile(grit_enroll,low) # lower limit of confidence interval
UL = quantile(grit_enroll,upp) # upper quantile for confidence interval
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
grit.enroll.fig <- hist(grit_enroll,breaks = 'FD',col = 'skyblue', 
                      xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                      main = '')

# Peer affiliation on currently enrolled
affiliation_enroll = (a2vec*enroll.b2vec)

LL = quantile(affiliation_enroll,low)
UL = quantile(affiliation_enroll,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
affiliation.enroll.fig <- hist(affiliation_enroll,breaks = 'FD',col = 'skyblue',
                             xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                             main = '')

# Self-control on currently enrolled
selfcontrol_enroll = (a3vec*enroll.b3vec)

LL = quantile(selfcontrol_enroll,low)
UL = quantile(selfcontrol_enroll,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
selfcontrol.enroll.fig <- hist(selfcontrol_enroll,breaks = 'FD',col = 'skyblue',
                             xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                             main = '')

# behavior on currently enrolled
behavior_enroll = (a4vec*enroll.b4vec)

LL = quantile(behavior_enroll,low)
UL = quantile(behavior_enroll,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
behavior.enroll.fig <- hist(behavior_enroll,breaks = 'FD',col = 'skyblue',
                          xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                          main = '')

# Spiritual on currently enrolled
spiritual_enroll = (a5vec*enroll.b5vec)

LL = quantile(spiritual_enroll,low)
UL = quantile(spiritual_enroll,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
spiritual.enroll.fig <- hist(spiritual_enroll,breaks = 'FD',col = 'skyblue',
                           xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                           main = '')

# Self identity on currently enrolled
identity_enroll = (a6vec*enroll.b6vec)

LL = quantile(identity_enroll,low)
UL = quantile(identity_enroll,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)
identity.enroll.fig <- hist(identity_enroll,breaks = 'FD',col = 'skyblue',
                          xlab = paste(conf, '% Confidence Interval', ' LL', LL4, ' UL', UL4),
                          main = '')

